Subspace-based line detection

ABSTRACT

A new signal processing method solves the problem of fitting multiple lines in a two-dimensional image. The Subspace-Based Line Detection (SLIDE) algorithm formulates the multi-line fitting problem in a special parameter estimation framework such that a signal structure similar to the sensor array processing signal representation is obtained. Any spectral estimation method can then be exploited to obtain estimates of the line parameters. In particular, subspace-based algorithms of sensor array processing (e.g., the ESPRIT technique) can be used to produce closed-form and high resolution estimates for line parameters. The signal representation employed in this formulation can be generalized to handle both problems of line fitting (in which a set of binary-valued discrete pixels is given) and of straight edge detection (in which one starts with a grey-scale image).

The invention relates generally to signal procssing for line detection,and more particularly the invention relates to fitting of multiple linesin a two-dimensional image.

BACKGROUND OF THE INVENTION

The background of the invention is described with reference to prior artpublications referenced in Appendix A attached hereto.

The problem of fitting lines to a given set of data points is awell-known scientific and engineering task. In computer vision, forexample, a digitized image may contain a number of discrete, `1` pixelslying on a `0` background, and the goal is the detection and parameterextraction of a number of straight lines that fit groups of collinear`1` pixels. A number of standard methods exist for such problems.Ordinary least squares methods, for example, seek to minimize the sum ofsquared vertical (or horizontal) distances of all points to the desiredline. However, the resulting solutions are coordinate-system dependent;moreover in ordinary least squares, all measurement error is attributedto only one coordinate of the data points. To avoid these problems, thetotal least squares method (TLS) [1, 2] uses a criterion in which thesum of squared normal distances of all data points to the line isminimized. However, there are two major limitations that are inheritedin all least squares approaches. One is the sensitivity of these methodsto outliers. In other words, although the least squares methods canefficiently handle measurement errors that appear as small deviations inthe point set, they cannot tolerate outliers; the existence of even arelatively small number of outliers tends to produce large residuals,leading to inadequate and unacceptable line fitting. The second mainlimitation of least squares methods lies in their inability to handlecases where there is more than one underlying line, especially crossinglines, in an image. The least squares method will result in an obviouslywrong answer in these cases; see FIG. 1 for a simple example.

Another now classical method of finding line parameters to fit a givendata point set, which is capable of handling multi-line problems, is theHough transform [3, 4, 5], which is a special case of the Radontransform [6]. In this approach, first a special transformation isapplied to all the points and then a two dimensional search isaccomplished to find the maxima in the transform plane. The lineparameters are ρ, the normal distance of the desired line from theorigin, and θ, the angle that the normal to the line makes with positivex-axis; ρ and θ are also the coordinates in the Hough transform outputplane. Each data point (x,y) is mapped under the Hough transform to thecurve

    ρ=x cos θ+y sin θ

in the ρ-θ plane. This equation also represents the line in the x-yplane that has a distance ρ to the origin and the normal to it makes anangle θ with the x axis. Therefore, all the points in the x-y planelocated on the line ρ₀ =x cos θ₀ +y sin θ₀ are mapped to curves in theρ-θ plane that all pass through the point (ρ₀,θ₀).

To fit a straight line to a set of data points, both ρ and θ axes haveto be quantized and hence a two-dimensional accumulator array beconstructed in the ρ-θ plane. The Hough transform equation is applied toeach point in the data set and the contents of all the cells in thetransform plane that the corresponding curve passes through areincreased by one. This is done until all the data points are transformedand then, an exhaustive search will be needed to locate a number ofmaxima points in the ρ-θ plane.

The Hough transform method is capable of handling a fairly high amountof noise contamination; however, it has certain drawbacks that maydrastically limit its use. One problem, mentioned earlier, is thequantization of the ρ and θ axes: to obtain an acceptable resolution,these coordinates need to be quantized finely, but reducing thequantization steps will broaden the peaks in the transform plane; inother words, since the actual measurement points for a line are notexactly collinear, the corresponding curves in the ρ-θ plane do not passthrough the same point. Certain procedures have been investigated tocompensate for this difficulty in the quantization step of the Houghtransform [7, 8, 9]. Secondly the burden of the two dimensionalexhaustive search in the ρ-θ plane for finding the maxima points is aserious drawback of the Hough transform, and this becomes more severewhen the number of lines to be fitted to a given data set is unknown. Insuch cases, the existence of local maxima may cause the search toprovide incorrect answers.

The present invention is directed to a new systematic multi-line fittingand straight edge detection procedure.

SUMMARY OF THE INVENTION

The Subspace-Based Line Detection (SLIDE) algorithm in accordance withthe invention, reformulates the line parameter estimation problem into asubspace fitting framework by applying a propagation (phase-preservingprojection) scheme on the image pixels. More specifically, SLIDE isbased on the introduction of a perfect mathematical analogy between theproblem of estimating angles of multiple lines in an image and theproblem of determining the directions of arrival of planarelectromagnetic waves impinging on an antenna array. Another analogy ismade between the problem of estimating the offsets of the lines and theproblem of time series harmonic retrieval. The direction-of-arrivalestimation and the harmonic retrieval problems have been extensivelystudied, e.g. for radar and sonar applications and time series analysis,and in particular, in the last decade several so called super-resolutionalgorithms (see [10],[11],[12],[13]) have been developed. In the linefitting problem, it is shown that the conditions for using a certaincomputationally efficient (ESPRIT) algorithm are met. We will introducethese analogies by adopting two signal generation scenarios from theimage pixels that we call the constant-μ propagation for line angleestimation, and the variable-τ propagation for line offset estimation.Another particularly appealing benefit of the analogy is that anefficient statistical technique from the field of sensor arrayprocessing can be used to provide an objective estimate of the number oflines in an image.

SLIDE yields closed-form and high resolution estimates of the lineparameters, as opposed to the standard method of Hough transform inwhich the line parameters are subject to limited resolution due to thechosen bin sizes and are found by implementing a search procedure.Computational complexity of SLIDE is an order of magnitude less thanthat of the Hough transform method. Also SLIDE does not require a largememory space, whereas the Hough transform method needs a large space tostore the output grid values.

In the new approach, estimation of line parameters is accomplished intwo phases. First, the directions (or angles) of the lines are estimatedand then estimates of the offsets are obtained. There exist importantapplications where only the line orientations are of concern. Forexample, when an autonomous vehicle is programed to follow a road, thescene that its camera pictures includes primarily (after preprocessing,)some `1` pixels representing the two sides of the road, as in FIG. 2(a),for example. By determining the orientations of these two lines, thevehicle can decide on its next movement. As another application, one canmention the task of precise rotational alignment of wafer and mask inmicrolithography for VLSI fabrication. In fact, it was the latterproblem that inspired the present work. An example of a preprocessedpicture from lithography is shown in FIG. 2(b).

We can now introduce a new formalism for the multiple line angleestimation problem that casts this problem into the general framework ofarray processing problems.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates an example of line fitting where the least squaresmethod results in a wrong answer.

FIGS. 2(a) and 2(b) illustrate two examples of practical use in whichonly the orientation of the lines is of concern; FIG. 2(a) is a typicalscene pictured by an autonomous vehicle commanded to move along a roadwhere by estimating the direction of the road sides, the vehicle decideson its next movement; and FIG. 2(b) illustrates, in precision wafer-maskrotational alignment, the direction of the feature-lines on the wafer isto be estimated.

FIG. 3 illustrates passive sensor array geometry.

FIG. 4(a) illustrates the image matrix and the hypothetical sensorslocated along the y-axis, and FIG. 4(b) illustrates measured signal atthe l^(th) receiver due to a single line.

FIG. 5 illustrates an overlapped subarray configuration in a uniformlinear array, the subarray separation Δ, equal to 1 times the elementspacing.

FIGS. 6(a)-6(c) are illustrations of measurement vectors in the lineangle and offset estimation problems: (a) image matrix, (b) projectionvector obtained in the constant-μ propagation for angle estimation, (c)projection vector obtained in variable-r propagation for offsetestimation, (top) before, and (bottom) after dechirping. Estimationresults: θ=-39.9905°, x0=100.9645 pixels.

FIGS. 7(a)-7(d) illustrate (a) image matrix, (b) the minimum descriptionlength function, (c) estimated line angles via different array procssingtechniques, and (d) regenerated image using the estimated lineparameters by SLIDE.

FIGS. 8(a) and 8(b) illustrate another example of line fitting by theproposed method, (a) original image with parallel and crossing lines,and (b) reconstructed image using the estimated line parameters:θ=30.0636°, -60,0041°, x0=98.5226, 199.6088, 300.1067, 400,0804,500.7668, and 100.7426 pixels.

FIGS. 9(a)-9(f) illustrate (a) an example of a noisy image, (b) Houghtransform output plane, (c)-(f) comparison of the results obtained bySLIDE and the Hough transform method, with dotted lines the results ofSLIDE and solid curves are cross sections of the Hough transform planewhere maxima occur, (c), (d) angle estimates, (e), (f) offset estimates.

FIGS. 10(a)-10(d) illustrate (a) original image showing part of asemiconductor wafer, (b) edge-enhanced image, and estimates of the linesangle (c), and offsets (d) by SLIDE and the Hough transform method.

FIGS. 11(a)-11(c) illustrate (a) aerial image, (b) high-contrast image,(c) detected straight roads superimposed on the contour plot of theaerial image. The angle estimates are -21.2685° and 24.7177° relative tothe vertical axis. The offset estimates are 48.5027 and 379.6986 pixels,respectively.

FIGS. 12(a)-12(c) illustrate an example of detecting linearly alignedobjects: (a) aerial airport image, (b) high-contrast image, (c)estimated support line of the airplanes superimposed on the contour plotof the aerial airport image. The angle estimate is 12.6081° relative tothe vertical axis. The offset estimate is 123.8449 pixels.

FIG. 13 illustrates corners of the rectangle which are detected byfitting lines to its sides.

DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS Estimation of LineAngles

SLIDE reformulates the multiple line angle estimation into the frameworkof sensor array processing by the virtue of introducing a propagationscenario for the image under consideration. In the so-called constant-μpropagation scenario [14], each pixel of the image in a row contributesto a received signal at a hypothetical sensor located in front of thecorresponding row of the image matrix. This contribution is in effectequivalent to a phase delay if the pixel is regarded to be a source ofelectromagnetic wave. More formally, let's assume that all pixels in theimage start to propagate narrowband electromagnetic waves with zeroinitial phases. Furthermore, assume that the waves estimated from pixelsin a given row of the image matrix are confined to travel only alongthat row towards the corresponding sensor. In this propagationenvironment, each straight line in the image will be in effectequivalent to a wavefront of a travelling planar wave, because the wavesof pixels on a line retain constant relative phase along paths parallelto that line. As we will see later, this propagation scenario createsmeasurements at the sensors similar to the measurements obtained in realantenna array processing.

We assume that a uniform linear array of sensors is located along thevertical axis of the image matrix (4(a)), each sensor receiving signalsonly from its corresponding row in the matrix according to the followingrule. If there are p nonzero pixels on the l-th row of the image matrixlocated on columns q₁, . . . q_(p) respectively, then the sensor infront of the l-th row will receive the signal ##EQU1## where μ is aconstant parameter to be chosen at our convenience [15]. Now assume thatthere is one straight line in an N×N image with offset x₀ and angle θ,as in FIG. 4. Provided that the line width is such that the line givesrise to only one nonzero pixel per row, and that there are no noisepixels, the measured signal at the l-th receiver will be

    z.sub.1 =e.sup.-jμx.sbsp.0 ·e.sup.jμl tan θ(2)

In a more general case, where there are d straight lines and also noisepixels in the image, the signal will be ##EQU2## where n_(l) includesthe effect of noise in the l-th row, which may consist of individualnoise pixels and also takes care of a possible displacement of a linepixel in that row.

Now let us define ##EQU3##

We will call this signal generation scenario the constant-μ propagation.This scenario is effectively encoding the line angles as frequencies ofcisoidal components in the measurement vector z=[z₀, . . . , z_(N-1)]^(T).

As Eq. (3) shows, the measurement vector z includes d cisoidalcomponents contaminated with noise. Any spectral estimation technique(including FET) can be applied in this stage to extract the frequenciesof the cisoids. However, since high resolution estimates of the anglesare desired, we pursue with exploiting the subspace-based techniques ofarray processing framework.

Using a window of length M, we can arrange the measurements in a matrixform as ##EQU4## where P=N+1-M and the new snapshot vector z_(p) εC^(M)×1 can be written as

    z.sub.p =A(θ)s.sub.p +n.sub.p                        (7)

In the above, A(θ)=[a(θ₁), . . . , a(θ_(d))] is an M×d matrix with k-thcolumn a(θ_(k))=[1.λ_(k), . . . , λ_(k) ^(M-1) ]^(T) with λ_(k) =e^(j)μtan θ.sbsp.k. This rearrangement of the data is called spatialsmoothing.

Combining Eqs. (5) and (6), we can write ##EQU5## contains the offsetsof the lines, and Φ=diag[λ₁, . . . , λ_(d) ].

Now we are ready to apply the well-studied array processing techniquesto the multi-line angle estimation problem. In sensor array processingframework, the goal is to estimate the directions of arrival θ of dimpinging planar wavefronts on the array of M sensors, given P snapshotsz_(i), i=1, . . . , P. There are many algorithms for solving thisproblem (see e.g. [10]-[13]), but for the present application aparticular one called TLS-ESPRIT appears to be the most appropriate.

In the ideal case of P→∞, the covariance matrix of the measurements hasthe following structure ##EQU6## and it is assumed that ##EQU7## Now letthe eigendecomposition of R_(zz) be ##EQU8## where λ₁ ≧λ₂ ≧ . . .≧λ_(M). The matrix E_(s) contains the d eigenvectors of R_(zz)corresponding to its d largest eigenvalues, i.e. E_(s) =[e₁, . . . ,e_(d) ]. The range space of E_(s) is called the signal subspace. Itsorthogonal complement is called the noise subspace and is spanned by thecolumns of E_(n) =[e_(d+1), . . . , e_(M) ]. It is easily seen from (11)that the smallest eigenvalue of R_(zz) has multiplicity M-d and is equalto the noise variance. Thus, Λ_(n) =σ² I. Moreover, it can be seen thatthe signal subspace is equal to span {A(θ)}, i.e. span {E_(s) }=span{A(θ)}, thus explaining the name. This is the basic observation in allsubspace-based array processing techniques, including ESPRIT.

In practice, only the same covariance matrix R_(zz).sbsb.P is available,and its eigendecomposition is ##EQU9## where now Λ_(n) ≠σ² I, exceptasymptotically as P→∞.

Determining the value of d is crucial here because by knowing d we canidentify the eigenvectors corresponding to the d largest eigenvalues ofR_(zz) as a basis for the estimate of the signal subspace. Under thestochastic signal assumption, the minimum description length (MDL)criterion is defined by ##EQU10## where λ_(i) are the eigenvalues ofR_(zz). The use of this criterion was first proposed in [16]. Theestimate of the number of sources d is given by the value of k for whichthe MDL function is minimized.

Now define E₁ as the submatrix E_(s) formed by rows 1 to M-1, and E₂ asthe submatrix of rows 2 to M. This is a special case of partitioning thereceiver array into two maximally overlapping subarrays with adisplacement Δ; here we have chosen Δ=1 (see FIG. 5). Next form theeigendecomposition of the 2d×2d matrix ##EQU11## Partition F into d×dsubmatrices, ##EQU12## and then let λ_(k), k=1, . . . , d be theeigenvalues of -F₁₂ F₂₂ ⁻¹. Finally the line angles are obtained as##EQU13##

Estimation of Line Offsets

In this section, we briefly discuss the method used in SLIDE algorithmfor estimating the line offsets.

After the angles have been estimated by the method of Sec. 3.1, theoffsets can be estimated by either of the three methods described inthis section.

As was mentioned above, the vector ##EQU14## contains the offsets of thelines. From the analogy with sensor array processing framework, solvingthe so-called signal copy problem yields estimates of values of elementsof s. This is done by solving the following problem

    s=argmin∥z-A(θ)s∥.sup.2            (20)

In a noise-free case, this criterion tries to find appropriate phasesfor cisoidal signals in order to fit their superposition to the givenmeasurement vector y [17]. It is easy to show that the distribution ofthe measurement noise is gaussian if the outliers are uniformlydistributed over the image [15]. Therefore, the above least squaresproblem provides the maximum likelihood estimate for the signal vector.Then, the offsets are computed via the relation ##EQU15##

However, the least squares technique is highly sensitive to outliers andits result rapidly degrades if the number of outliers increases. Also inmany practical situations, the image may contain only a part of astraight line. In such cases, the above least squares estimate will bebiased and unsatisfactory. There also exists a degenerate case in usingthis technique for offset estimation. Namely, since each straight lineis assumed to represent a travelling planar wavefront, it can be shownthat the phases of the signals originating from parallel lines cannot beresolved at the sensors. In other words, the received signal will be acoherent superposition of the two or more parallel wavefronts.

Another method for estimating the offsets is to project the image alongeach of the detected angles and search for a peak in the resultingprojection. This method is robust to outliers but suffers from limitedresolution problem as we have to choose bin sizes for the projection.

In order to obtain closed-form and high resolution estimates for theoffsets, SLIDE proposes another propagation scenario that leads toencoding of the line offsets as frequencies of cisoidal components inthe measurement vector. We will call this propagation scheme thevariable-τ propagation for reasons that will become clear shortly. Thepropagation part of this procedure can be done in parallel with theconstant-μ propagation that encodes the angles, but the offsetestimation procedure that follows the new propagation scenario has to bedone after the line angles have been estimated.

Assume that there are d₁ parallel lines in the image with angle θ₁.Performing the same propagation scheme of above to the pixels on thel-th row and replacing the name of propagation parameter μ with τ, thesignal received at the l-th sensor is ##EQU16## Now we let τ vary alongthe rows of the matrix. For simplicity of the analysis, here we choose τto be proportional to the row number l

    τ=l/α                                            (23)

where α is a constant. Thus, the measurement vector z will now haveelements that are functions of τ ##EQU17## As we observe in Eq. (24), inthe new measurements τ appears in quadratic form with the angle and inlinear form with the offsets. If there is only one line in the image (d₁=1), then the sequence z.sub.τ will be of the form

    z.sub.τ =e.sup.jαr.spsp.2.sup.tan θ.sbsp.1 ·e.sup.-jτx.sbsp.0 +n.sub.τ, τ=0, . . . , α(N-1)                                              (26)

which has the form of a chirp signal. In general, z.sub.τ in Eq. (24)represents the sum of d₁ chirp signals with the same quadratic term butwith different linear terms. The quadratic term can be factored out ofthe summation and as before be called a.sub.τ (θ₁) as in Eq. (25). Sincewe know the angle θ₁ from the previous analysis (constant-μpropagation), we can divide z.sub.τ by a.sub.τ (θ₁) and get anothersequence w.sub.τ which is free from the quadratic term in τ. This isbasically equivalent to dechirping the measurement sequence. The newmeasurement vector w has elements of the following form ##EQU18## wheren'.sub.τ =n.sub.τ /a.sub.τ (θ₁) is the new additive noise term.

The sequence w.sub.τ in Eq. (28) is a superposition d₁ cisoids plusnoise. The frequency of each of the cisoids is proportional to one ofthe line offsets. Therefore, we have reformulated the offset estimationproblem as a time series harmonic retrieval problem. By exploiting theeigenstructure of the covariance matrix of the measurements and subspacefitting techniques of sensor array processing, high resolution estimatesof the offsets can be obtained. Also the number of parallel lines foreach angle can be estimated using the minimum description lengthcriterion.

In general, an image may contain parallel and non-parallel lines. Inorder to analyze the nature of the measurement vector w in that case,let's assume that there are d₁ parallel lines in the image with angle θ₁and d₂ parallel lines with angle θ₂. We order the offsets of all thelines such that ##EQU19## correspond to the lines with angle θ₁ and##EQU20## correspond to the lines with angle θ₂. The z vector will haveelements ##EQU21## Now consider the estimation of the offsets of the newlines having angle θ₁. As before, by dividing the sequence z.sub.τ bya.sub.τ (θ₁), we dechirp the components of the measurement sequence thatare contributions from angle θ₁. This results in ##EQU22## The termscorresponding to angle θ₁ have been dechirped and transformed to purecisoids as before. The terms corresponding to angle θ₂ however, havebeen transformed to another chirp form. The basic assumption forapplying the ESPRIT algorithm to the harmonic retrieval problem is thatthe components s(t) of the underlying analog signal be narrowband, i.e.for each component we can write s(t-t₀)≈e^(-j)ωt.sbsp.0 s(t). If onlythe pure cisoids corresponding to the liens with angle θ₁ were present,this condition would be directly satisfied. The sequencew.sub.τ.sup.(θ.sbsp.1.sup.) in Eq. (33) contains both pure cisoids andchirp components. However, the energy of the chirp components isdistributed over the frequency axis and they can be considered as noisein the analysis. Therefore, the harmonic retrieval procedure will beable to detect the pure cisoids corresponding to the angle underconsideration. More clarification and examples for this point will bepresented later. In a general case, the variable-τ propagation needs tobe performed only once, and the dechirping and harmonic retrievalprocedures have to be done for each of the detected angles.

The procedure of above is then used to rearrange the sequence w.sub.τinto a Hankel matrix W_(P) as in Eq. (6) and then to implement thesubspace-based estimation techniques such as ESPRIT to solve for thefrequency components in w.sub.τ which correspond to the line offsets.

Practical Issues

Generalizing to Grey-Scale Images

The formalism introduced above can also be generalized to handle theproblem of straight edge detection in grey-scale images. Namely, we canstart from a grey-scale image and detect parameters of straight lines inthat image. This is accomplished by simply assigning an amplitude to thepropagated signal from each pixel proportional to the grey-scale valueof that pixel. A primitive edge enhancement procedure may be used firstin this case to attenuate background contributions. This edgeenhancement can in general be as simple as a shirt-and-subtract or acontrast enhancement procedure applied to the image.

Mathematically, if the edge strength in the first row (after edgeenhancement) is represented by a function f(x)=g(x-x₀), then the sensorin front of the first row will receive ##EQU23## where G(μ) is equal tothe Fourier transform of g(x) at a single frequency point, and is aconstant. In binary images we have G(μ)=1. The edge strength in row l(in a one-line and noise-free image) will be equal to g(x-x₀ +l tan θ),and the receiver in front of that row will read ##EQU24## The onlydifference between the above equation and Eq. (2) is the appearance of aconstant (G(μ)) as the amplitude. This treatment shows that identicalprocedures as those that are proposed for line fitting in binary imagescan be performed on edge-enhanced grey-scale images.

A case of interest in line fitting in grey-scale images isnon-uniformity in edge strength over different image lines.Mathematically, the deviation in value from the nominal (or average)edge strength can be lumped into the noise term for each line.Simulation results show that as long as the signal to noise ratioremains in a reasonable level, the algorithm is still able to detect theedge in these cases.

Aliasing Issues

As was mentioned above, the solution is provided in two steps. First,the line directions are estimated and then an estimate of the offsets isobtained. We start with the estimation of line angles θ_(k), k=1, . . .d. Using the TLS-ESPRIT algorithm, the angles are estimated via therelation ##EQU25## where Δ is the displacement between the two subarraysused in the ESPRIT algorithm. In the above equation, normalization isdone in order to compensate for possible perturbations in magnitude ofλ_(k) =e^(j)μΔ tan θ.sbsp.k from unity. The range of non-ambiguousangles is determined by the following formula ##EQU26## Values of μ andΔ can be chosen appropriately in order to avoid the ambiguity problem.It is also worth mentioning that since the value of μ is not dictated bythe physics as it is the case in sensor array processing, improvedresolution can be achieved by increasing its value. To be more specific,from Eq. (37) we must choose ##EQU27## in order to avoid ambiguity. Inpractice however, due to the fact that the image is in fact discretized(sampled), the value of μ cannot be chosen higher than π. For angleestimation stage, μ is better to be chosen to have a value around unityto provide both non-ambiguity and good resolution. In estimating theoffsets, the value of μ should be chosen such that the line offsets donot become ambiguous. In other words, since the received signal due to apixel at a receiver has a period equal to 2π, and offsets can havevalues in the range [0, N-1], μ has to be chosen such that |μ|≦2π/N.Therefore, it is recommenced that the propagations performed in angleand offset estimation steps be performed with different values of μ.

The aliasing issue arises also in the line offset estimation problem.Here, the factor that needs to be controlled is the step size chosen forthe sequence τ, which we call δτ. Larger values of δτ result in moreperiods of the cisoids in the data records which is desirable ifaliasing does not occur. In general, the value of δτ should satisfy##EQU28## in order that aliasing does not happen. However, if anestimate of the maximum possible value for all line offsets isavailable, the N in the above formula can be replaced by that value andδτ can be increased accordingly.

Combining Different Estimation Results

Another flexibility of the proposed algorithm that can be exploited toimprove the signal to noise ratio of the angle estimation procedure isthe possibility of doing estimation using hypothetical sensors along thex-axis as well and then combining the results of the two estimations.

In offset estimation procedure, the mentioned combination of estimatesalong the x-axis and the y-axis can help to improve the results.Moreover, we can observe that in this case the lines that have largeroffsets are transformed to cisoids with higher frequencies. Thesefrequencies can normally be detected more accurately due to the presenceof many cycles of the signal in the data record. Therefore, thepropagation can be performed along both positive and negative directionsof the coordinate axes, and in each case the more accurate results bechosen according to the location of the estimated offsets. For example,with propagating along negative x-axis, offset estimates with valueslarger than N/2 are more accurate, while with propagating along positivex-axis, the estimates that are less than N/2 are preferred.

Partial Extent Lines

Another practical issue arises when there are line segments in the imagerather than long lines. The spatial smoothing procedure provides a levelof robustness for the process of angle estimation by sliding a smallsize window over the measurement vector. This way, a line with a supportlong enough to be included in many windows will be adequatelyrepresented in the covariance matrix. However, small support linescannot be represented adequately and won't be resolved in theestimation. A remedy for such a case is to divide the image into severalsmaller sub-images and perform the analysis on each sub-image. Theresults of individual sub-image estimations can be put together toobtain parameters of both long and short lines in the image. In order toclarify why lines that have partial support in the image can bedetected, the following observation is also worth mentioning. The basicassumption for applying the subspace techniques to a set of subsequentsnapshots is that there exists a certain phase relation betweensubsequent snapshots. The matrix Φ in the ESPRIT algorithm is in fact adiagonal matrix containing these phase relations for different signalcomponents. In the SLIDE algorithm, this phase relation is retainedbetween subsequent subvectors made from measurements that result fromlines with partial extent. More specifically, snapshots (subvectors)that are fully covered by the lateral extent of the line clearly satisfythe phase relationship. Snapshots that come from the outside of the linesupport have zero contribution to the estimate (although in the presenceof noise they tend to lower the SNR). Finally, for snapshots that arepartially covered by the line lateral extent, the phase relation ispreserved between each two subsequent snapshot in all but two pointsthat correspond to the ends of the line support. This is the case simplybecause the same phase relationship that exists for the nonzero part ofthe snapshots can be assumed to exist for the identically-zero parts ofthem. Another viewpoint about this observation follows shortly.

Choice of the Window Size

In the procedure of windowing the long N×1 measurement vector z (and w)to produce the matrix Z_(P) (and W_(P)), the values of P and M need tobe chosen such that P plays the role of number of snapshots and M is thedimension of the snapshots. Obviously, for the sample covariance matrixto approximate the true covariance matrix, we need to choose a large P.On the other hand, choosing the dimension M of the new vectors largehelps to smooth out the eigenvalues corresponding to the noise subspaceof the covariance matrix. Therefore, in order to make the best use ofthe measurement vector, sub-vectors are chosen in such a way to maintainmaximum overlap with their adjacent neighbors. With this choice, therewill be a fixed relation between the number of snapshots P, and thedimension of the new measurement vector M, namely P+M=N+1. It shouldalso be mentioned that choosing the dimension of new vectors to be verylarge increases the dimension of the sample covariance matrix. Thisresults in dramatic increase in computation load if one tries to obtainthe whole eigendecomposition of the sample covariance matrix; however,if the so called Fast Subspace Decomposition (FSD) technique of Xu andKailath [18] which is based on the Lanczos algorithm [19, 1] is used toextract only the signal subspace of the sample covariance matrix, thecomputation load will increase linearly with M.

Computational Complexity of Slide vs the Hough Transform

In order to compare the computational effort needed in SLIDE with thatof the Hough transform method, two cases are considered here. In thefirst case, it is assumed that the image is binarized and that there isone pixel in each row of the N×N image matrix for every one of the dlines. Without the presence of any outliers, the computational effortneeded for applying the Hough transformation to all pixels will be ofO(N² d) provided that we choose to discretize the θ axis to N bins; insome applications, we may have prior information about the possiblerange of the angles, and discretize only that range into N bins. Inaddition to the above, a search procedure is needed for estimating thelocation of peaks in the transform plane. On the other hand, SLIDE needsO(Nd) calculations for the propagations and O(N M²) for the completeeigendecomposition of the sample covariance matrix. However, if thestructure of the data matrix is exploited and also only the signalsubspace of the sample covariance matrix is extracted as in the newmethod of Xu and Kailath [18] (note that this subspace contains allneeded information for estimating the parameters), then the totalcomputation load of SLIDE will be further decreased to O(N(M+2d)). Inpractice, the value of M is chosen around √N, so the computationalcomplexity of SLIDE in this case is O(N^(3/2)) compared to O(N² d) forthe Hough transform method. In the second case, it is assumed that allthe N² pixels in the image participate in the computations (e.g. we havea grey-scale image). In this case, the Hough transform method requiresO(N³) calculations plus search. The computational effort for SLIDE inthis case will be O(N²) for the propagations and again O(N(M+d)) for therest of the calculations. So the computational complexity of SLIDE inthis case is O(N²) as compared to O(N³) for the Hough transform method.This reduction in computational complexity could allow real-timeimplementation of SLIDE.

Summary of the Slide Algorithm

In this section, an outline of the SLIDE algorithm as applied forestimating angles and offsets of straight lines and edges in an image ispresented. Again, the image size is assumed to be N×N.

If the image is binary, go to the next step. If the image is grey-scale,produce an edge-enhanced image from the digitized grey-scale image. Thiscan be done by a simple shift an subtract (one-pixel gradient) operatoror a contrast enhancement procedure.

Perform the constant-μ propagation as discussed above on the image andobtain the measurement vector z. In edge-enhanced grey-scale images,this can be done by multiplying each row of the image matrix by thevector [1,e^(-j)μ, e^(-j2)μ, . . . , e^(-j)(N-1)μ ]^(T). However, inbinary images savings in computations can be made by performingcalculations only on nonzero pixels of the image.

In this stage, any spectral estimation method can be applied to estimatethe cisoidal components of z. However, in order to obtain highresolution and closed-form estimates, use subspace-based techniques ofarray processing as follows.

Rearrange the measurement vector using Eq. (6) such that P snapshotvectors of size M are produced.

Compute the sample covariance matrix of the snapshots as in Eq. (15);and perform its eigendecomposition. By applying the MDL method, thenumber of present angles can be estimated.

Estimate the angle by any subspace-based method. An outline of how toexploit the ESPRIT algorithm is presented above. ESPRIT provides highresolution estimates without the need for a search procedure in theparameter space.

After knowing the angles of the lines, estimate the offsets by any ofthe three methods that follow. Using least-squares minimization as inEq. (20), or projecting the image along each detected angle andsearching for peaks, or applying the variable-τ propagation stage ofSLIDE as described in Sec. 3.2 and summarized below.

For exploiting the offset estimation part of SLIDE, perform thevariable-τ propagation (as described in Sec. 3.2), by multiplying eachline of the image matrix by a vector as follows. Numbering the rows ofthe image matrix from 0 to N-1, multiply the l-th row by the vector [1,e^(-jl/)α, e^(-j2l/)α, . . . , e^(-j)(N-1)l/α ]^(T) to obtainmeasurement vector z. (Again, in binary images computational savings canbe made by only performing calculations only on nonzero pixels.)

If there are d different angles in the image, perform the followingstages for dechirped measurements for each angle obtained from themeasurement vector.

Divide each element of the measurement vector by the correspondingelement response value as in Eq. (28) to get the dechirped measurementvector w.

Rearrange the elements of w the same way the measurement vector z wasrearranged in angle estimation stage to get P snapshots each of lengthM.

Compute the sample covariance matrix of the new snapshots and performits eigendecomposition. Use the a priori knowledge or MDL to decide onthe number of present parallel lines.

Exploit the ESPRIT algorithm (or other spectral estimation methods) toget estimates of the offsets.

Experimental Results

In this section, several experimental results of applying the SLIDEalgorithm to estimating the parameters of lines in images are provided.In order to make comparisons, the method of Hough transform which isessentially a projection technique is also applied to images and theresults of both approaches are presented. It is worth re-emphasizinghere that the estimates provided by the SLIDE algorithm are highresolution while the Hough transform (and projection) results arelimited in resolution to the bin sizes chosen for the parameters. Also,SLIDE is computationally much more efficient than the Hough transformand there is no need for a search procedure in the parameter space orcurve reading in its implementation.

FIG. 6 illustrates the measurement vectors obtained in constant-μpropagation and variable-τ propagation. FIG. 6(a) shows an imagecontaining a straight line and outliers. In FIG. 6(b) the measurementvector obtained in constant-μ propagation is shown. In FIG. 6(c), themeasurement vector obtained in variable-τ propagation is shown (top)before, and (bottom) after dechirping. The line angle and offsetestimation problems are transformed to the problem of estimatingfrequencies of cisoids contaminated in noise.

In FIG. 7(a), an image containing five nearly collinear pixel sets isshown. The line angles are {-40, -35, 25, 30, 50} degrees and the lineoffsets are [25, 45, 150, 175, 190} pixels, respectively. Afterperforming the constant-μ propagation and computing the samplecovariance matrix of the measurements, the minimum description lengthfunction is calculated. This function is shown in FIG. 7(b); thelocation of its minimum is an estimate of the number of lines in theimage. In FIG. 7(c), the results of applying different sensor arrayprocessing techniques to this example are presented; the ESPRITalgorithm requires no search procedure and is by far the fastest one touse. The estimated angles by SLIDE are {-40.0424, -35.0399, 25.2412,30.2075, 50.0380} degrees and the estimated offsets are {25.7130,44.9882, 151.1676, 174.0602, 190.0925} pixels, respectively. Using thesevalues, the image of FIG. 7(d), is regenerated.

Another example is presented in FIG. 8 in which parameters of theparallel and crossing lines in an image are estimated.

To compare the efficiency and speed of the proposed technique with thestandard method of Hough transform, we consider the rather noisy imageof FIG. 9(a). This image is 400×400 pixels and was generated to containtwo lines with random angles θ=19.8418°, 27.2582° and random offsets x₀=186.8439, 277.2624 pixels, respectively. Note that due todiscretization effects, the actual values may be slightly different fromthe above. For applying the Hough transform, we need to choose a searchrange for the angle and a bin size for discretizing the normal distance(ρ) of lines from the origin. For this example, a uniform range ofangles was chosen between -35° and 35° with 0.1° steps. The bin size forρ was chosen to be 0.5 pixels. Note that ρ is the normal distance oflines from the origin, and the corresponding x-axis offset is computedby ##EQU29## FIG. 9(b) shows the Hough transform output plane, where thecoordinates of the peaks yield estimates for θ and ρ. Cross sections ofthe Hough transform plane at the locations of the peaks are depicted inFIG. 9(c)-(f) together with estimates obtained by SLIDE. The numericalvalues are as follows. The estimates by SLIDE are θ=-19.7485°, 27.2491°and x₀ =187.1408, 277.6463 pixels. Performing a search in the Houghtransform output plane, yields the angle estimates as θ=-19.8°, 27.2°and the offset estimates as x₀ =187.6, 277.7 pixels. It is seen thatboth techniques yield rather accurate estimates, however, the effort ofcomputations of SLIDE is about an order of magnitude (400 times in thisexample) less than that of the Hough transform. Also, the Houghtransform method needs a large memory space to store and accumulate theoutput.

As is discussed above, a straight-forward generalization of the proposedalgorithm for being applied to grey-scale images can be readily derived.Consider the grey-scale image shown in FIG. 10(a). First, a primitiveedge enhancement was done on this image in order to attenuate thebackground contributions. This was done by shifting the image matrix byone pixel, subtracting the result from the original matrix, and takingthe absolute value of the resulting matrix. The image of FIG. 10(b)shows the result of this edge enhancement procedure. Then, SLIDE and theHough transform method were applied to this image. FIGS. 10(c), (d) showthe angle and offset estimation results, respectively.

In FIG. 11, the result of applying SLIDE to a real aerial image ispresented. FIG. 11(a) shows the original aerial image. The high-contrastimage of FIG. 11(b) is used for line detection. In FIG. 11(c), theresults of line detection are superimposed on the contour plot of theimage.

FIG. 12 is an example of detecting linearly aligned objects by SLIDE.FIGS. 12(a), (b) show the original and high-contrast aerial images inwhich several airplanes are aligned along a certain direction which wemight like to determine. This is an illustrative example of non-uniformedge strength over different image rows. Applying the SLIDE algorithm tothe high-contrast image results in line parameters that are used togenerate the dotted line in FIG. 12(c). As this figure suggests, thedirection along which the airplanes lie seems to be well estimated bythe algorithm.

In FIG. 13, SLIDE is applied to detection of the sides (and hence thecorners) of a rectangle. Although the sides of the rectangle have alimited extent in the image (and in the measurement vector), theestimates of the corners are reasonably accurate.

While the invention has been described with reference to specificembodiments, the description is illustrative of the invention and is notto be construed as limiting the invention. Various modifications andapplications may occur to those skilled in the art without departingfrom the true spirit and scope of the invention as defined by theappended claims.

APPENDIX A References

[1] G. H. Golub and C. F. Van Loan, Matrix Computations, Johns HopkinsUniversity Press, Baltimore, Md., 1984.

[2] S. Van Huffel and J. Vandewalle, "The Total Least Squares Technique:Computation, Properties and Applications", In E. F. Deprettere, editor,SVD and Signal processing: Algorithms, Applications and Architectures,pages 189-207. Elsevier Science Publishers, 1988.

[3] P. Hough, "Method and Means for Recognizing Complex Patterns", U.S.Pat. No. 3,069,654, 1962.

[4] R. O. Duda and P. E. Hart, "Use of the Hough Transform to DetectLines and Curves in Pictures", Communications of the ACM, 15:11-15,1972.

[5] A. K. Jain, Fundamentals of Digital Image Processing, Prentice-Hall,Englewood Cliffs, N.J., 1989.

[6] S. R. Deans, The Radon Transform and Some of Its Applications, JohnWiley & Sons, New York, N.Y., 1983.

[7] W. Niblack and D. Petkovic, "On Improving the Accuracy of the HoughTransform: Theory, Simulations and Experiments", In Proc. IEEE Conf. onComputer Vision and Pattern Recognition CVPR'88, pages 574-579, 1988.

[8] E. R. Davies, Machine Vision: Theory, Algorithms, Practicalities,Academic Press, London, 1990.

[9] N. Kiryati and A. M. Bruckstein, "On Navigation Between Friends andFoes", IEEE Trans. on Pattern Analysis and Machine Intelligence,13(6):602-606, 1991.

[10] R. O. Schmidt, "Multiple Emitter Location and Signal ParameterEstimation", In RADC Spectrum Estimation Workshop, Griffiss AFB, NY,1979.

[11] A. Paulraj, R. Roy, and T. Kailath, "Estimation of SignalParameters by Rotational Invariance Techniques (ESPRIT)", In Proc. of19th Asilomar Conference on Circuits, Systems and Comp., 1985.

[12] M. Viberg, B. Ottersten, and T. Kailath, "Detection and Estimationin Sensor Arrays Using Weighted Subspace Fitting", IEEE Trans. on SP,39(11):2436-2449, Nov. 1991.

[13] R. Roy and T. Kailath, "ESPRIT: Estimation of Signal Parameters viaRotational Invariance Techniques", IEEE Trans. on ASSP, 37(7):984-995,July 1989.

[14] H. K. Aghajan and T. Kailath, "SLIDE: Subspace-based LineDetection", In Proc. of IEEE ICASSP, Minneapolis, Minn., to appear,April 93.

[15] H. K. Aghajan and T. Kailath, "Sensor Array Processing Techniquesfor Super Resolution Multi-Line Fitting and Straight Edge Detection",IEEE Trans. on Image Processing, to appear, 93.

[16] M. Wax and T. Kailath, "Detection of Signals by InformationTheoretic Criteria", IEEE Trans. on ASSP, 33(2):387-392, April 1985.

[17] H. K. Aghajan and T. Kailath, "A Subspace Fitting Approach to SuperResolution Multi-Line Fitting and Straight Edge Detection", In Proc. ofIEEE ICASSP, pages III:121-124, San Francisco, Calif., 1992.

[18] G. Xu and T. Kailath, "A Fast Algorithm for Signal SubspaceDecomposition and Its Performance Analysis", In Proc. of IEEE ICASSP,pages 3069-3072, Toronto, Canada, 1991.

[19] B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice-Hall,Englewood Cliffs, N.J., 1980.

What is claimed is:
 1. A method of estimating line parameters includingangles and offsets of straight lines in a two-dimensional imagecomprising the steps ofa) producing an edge-enhanced image from agrey-scale image, b) performing a constant-μ propagation on the image toobtain a measurement vector z, c) estimating the cisoidal components ofsaid measurement vector z including the number of angles present, and d)estimating the offsets of lines for each angle.
 2. The method as definedby claim 1 wherein step c) includesrearranging the measurement vectorsuch that P snapshot vectors of size M are produced, computing the samecovariance matrix of said snapshots, and performing eigendecompositionof said sample covariance matrix to estimate the number of anglespresent.
 3. The method as defined by claim 1 wherein step d) includesusing a least-squares minimization.
 4. The method as defined by claim 1wherein step d) includes projecting images along each detected angle andsearching for peaks.
 5. The method as defined by claim 1 wherein step d)includes applying a variable τ propagation stage of a subspace-basedline detection algorithm.